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Abstract 

We  propose  a  new  formulation  using  the  closest  point  mapping  for  integrat¬ 
ing  over  smooth  curves  and  surfaces  with  boundaries  that  are  described  by  their 
closest  point  mappings.  Contrary  to  the  common  practice  with  level  set  meth¬ 
ods,  the  volume  integrals  derived  from  our  formulation  coincide  exactly  with  the 
surface  or  line  integrals  that  one  wish  to  compute.  We  study  various  aspects 
of  this  formulation  and  provide  a  geometric  interpretation  of  this  formulation 
in  terms  of  the  singular  values  of  the  Jacobian  matrix  of  the  closest  point  map¬ 
ping.  Additionally,  we  extend  the  formulation  -  initially  derived  to  integrate 
over  manifolds  of  codimension  one  -  to  include  integration  along  curves  in  three 
dimensions.  Some  numerical  examples  are  presented. 


1  Introduction 

This  paper  provides  simple  formulations  for  integrating  over  manifolds  of  codimensions 
one,  or  two  in  M3 .  when  the  manifolds  are  described  by  functions  that  map  points  in 
Mn  (n  =  2, 3)  to  their  closest  points  on  curves  or  surfaces  using  the  Euclidean  distance. 
The  idea  for  the  formulation  originated  in  [7]  where  the  authors  proposed  a  formulation 
for  computing  integrals  of  the  form 


f  v(x.(s))ds,  (1) 

JdQ 

in  the  level  set  framework,  namely  when  the  domain  0  is  represented  implicitly  by 
the  signed  distance  function  to  its  boundary  Typically  in  a  level  set  method 

[11,  12,  14],  to  evaluate  an  integral  of  the  form  of  (1)  where  dQ  is  the  zero  level  set 
of  a  continuous  function  (/?,  it  is  necessary  to  extend  the  function  v  defined  on  the 
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boundary  <9D  to  a  neighborhood  in  MV  The  extension  of  u,  denoted  h,  is  typically  a 
constant  extension  of  v.  The  integral  is  then  approximated  by  an  integral  involving  a 
regularized  Dirac-5  function  concentrated  on  <9D,  namely 

/  v(x.(s))ds  «  /  u(x)5e((/;(x))| V(/p(x)|dx. 

Jon  J Mn 

Various  numerical  approximations  of  this  delta  function  have  been  proposed,  see  e.g. 
[3,2,15,16,17], 

In  [7],  with  the  choice  of  <p  =  dgn  being  a  signed  distance  function  to  9Q,  the 
integral  (1)  is  expressed  as  an  average  of  integrals  over  nearby  level  sets  of  dgn,  where 
these  nearby  level  sets  continuously  sweep  a  thin  tubular  neighborhood  around  the 
boundary  9f2  of  radius  e.  Consequently,  (1)  is  equivalent  to  the  volume  integral  shown 
on  the  right  hand  side  below: 


/  v(x(s))ds  =  /  v(x*)J(x-1dga)S€(ddn(x))dx, 
JdQ  J Rn 


(2) 


where  Se  is  an  averaging  kernel,  x*  is  the  closest  point  on  <90  to  x  and  J(x;  dm) 
accounts  for  the  change  in  curvature  between  the  nearby  level  sets  and  the  zero  level 
set. 

Now  suppose  that  <90  is  a  smooth  hypersurface  in  M3  and  assume  that  x  is  suffi¬ 
ciently  close  to  0  so  that  the  closest  point  mapping 


x*  =  Pgn(x)  =  argminyeaa|x  -  y| 

is  continuously  differentiable.  Then  the  restriction  of  Pgn  to  <90);  is  a  diffeormorphism 
between  9 and  90,  where  <90^  :=  {x  :  dgn(x)  =  r/}.  As  a  result,  it  is  possible  to 
write  integrals  over  90  using  points  on  90y  as: 

f  v(x)dS  =  f  v(x*)J(x.-,r])dS, 

J dQ  J dQv 

where  J(x,rj )  comes  from  the  change  of  variable  defined  by  Pgn  restricted  on  90f/. 
Averaging  the  above  integrals  respectively  with  a  kernel,  Se,  compactly  supported  in 
[— e,  e],  we  obtain 


v(x)dS 


v(x*)J(x-,r))dS  dr). 


Formula  (2)  then  follows  from  the  coarea  formula  [5]  applied  to  the  integral  on  the 
right  hand  side. 

This  paper  is  motivated  by  the  recent  success  in  the  closest  point  methods  for 
solving  partial  differential  equations  on  surfaces  [8,  9,  10,  13]  as  well  as  the  need  to 
process  data  sets  that  contain  unstructured  points  sampled  from  some  underlying 
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surfaces.  In  the  following  section,  we  show  that  in  three  dimensions  the  Jacobian  J 
in  (2)  is  the  product  of  the  first  two  singular  values,  G\  and  <x2,  of  the  Jacobian  matrix 
of  the  closest  point  mapping  u1^;1 ;  namely, 

/  v(x(s))ds  =  /  u(P0o(x))Je(d9o(x))TT(Ti(x)dx. 

Jan  Jr3  -=1 

To  motivate  the  new  approach  using  singular  values,  we  consider  Cartesian  coordinate 
systems  with  the  origin  placed  on  points  sufficiently  close  to  the  surface,  and  the 
^  direction  normal  to  the  surface.  Thus  the  partial  derivatives  of  the  closest  point 
mapping  in  the  z  direction  will  yield  zero  and  the  partial  derivative  in  the  other  two 
directions  naturally  correspond  to  differentiation  in  the  tangential  directions.  Thus 
we  see  that  one  of  the  singular  values  should  be  0  while  the  other  two  are  related  to 
the  surface  area  element.  We  also  derive  a  similar  formula  for  integration  along  curves 
in  three  dimensions  (codimension  2).  The  advantages  of  this  new  formula  include  the 
ease  for  constructing  higher  order  approximations  of  J  via  e.g.  simple  differencing, 
even  in  neighborhoods  of  surface  boundaries  where  curvatures  become  unbounded. 

For  clarity  of  the  exposition  in  the  rest  of  the  paper,  we  will  now  denote  the  distance 
function  simply  by  d. 


2  Integration  using  the  closest  point  mapping 

In  this  section,  we  relate  the  Jacobian  J  in  (2)  to  the  singular  values  of  the  Jacobian 
matrix  of  the  closest  point  mapping  from  R2  or  R3  to  T,  where  T  denotes  the  curves 
or  surfaces  on  which  integrals  are  defined.  We  assume  that  in  three  dimensions,  if  T 
is  not  closed,  it  has  smooth  boundaries. 


2.1  Codimension  1 

We  consider  a  C2  compact  curve  or  surface  T  that  can  either  be  closed  or  not.  If  T 
is  closed,  then  it  is  the  boundary  of  a  domain  fl  so  that  T  can  be  denoted  (XI.  If  T 
is  not  closed,  we  assume  that  it  has  smooth  boundaries.  We  define  d:Rn4RU  {0} 
to  be  the  distance  function  to  T  and  Pp  to  be  the  closest  point  mapping  Pr  :  Rn  i->-  T 
(for  n  =  2,3)  defined  as 

|Pp(x)  —  x|  =  min  |y  —  x|.  (3) 

We  let  do  be  the  distance  function  to  T  if  it  is  open  and  ds  be  the  signed  distance 
function  to  T  =  cXl  if  it  is  closed.  The  signed  distance  function  is  defined  as 


ds(x) 


infyGftc  |x  -  y|  if  x  G  S2, 
-  infyGn  |x  -  y|  if  x  G  flc. 
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Then  we  define  d  as  follows: 


d(x) 


{do(x)  if  T  is  open, 
ds(x)  if  T  is  closed. 


(4) 


Lemma  1.  Let  d  be  the  distance  function  to  Y  defined  in  (4).  For  \rj\  sufficiently  close 
to  0,  the  Gaussian  curvature  at  a  point  on  the  p  level  set  Yv  :=  {f  :  d(f)  =  p]  can  be 
expressed  as 


Gn  dXXdyy  +  dXXdZZ  +  dyyd 


yy^zz 


d2 

11  xy 


d: 


2 

‘xz 


d2 

ayz- 


(5) 


Proof.  Starting  with  the  definition  of  the  Gaussian  curvature  G  for  a  surface  (see  [6]), 
we  can  obtain  an  expression  for  the  Gaussian  curvature  of  its  p- level  set  in  terms  of  d 
as 


G=  (Vd,adj(Hess(d))Vd) 

dx(dyydzz  dyf)  ~\~  dy(dxxdzz  dxz )  +  dz(dxxdyy  dXy ) 

4~  ‘2\dxdy(KdxzdyZ  dXydZz)  -\-  dydz(^dXydxz  dyZdXJfj 

T  dXdZ(dXydyZ  dXZdyyf\. 


(6) 


We  show  that  this  expression  is  the  same  as  (5)  by  rearranging  the  terms  above  and 
using  the  fact  that  close  to  Y  the  distance  function  satisfies  |Vd|  =  1.  First  we 
rearrange  the  terms  in  G: 


G  dxdyydzz  dydxxdzz  d zdxxdyy  dxdyZ  dydxz  dzdXy 

T2  \dxdy  ( dxzd  yz  dXydZZ^j  +  dydZ(dXydXZ  dyZdXfj  ~\~  dXdZ(dXydyZ  dXZdyyf  J 

and  rewrite  each  of  the  first  six  terms  in  terms  of  |Vd|2,  e.g. 

dxdyydZZ  dyydZZ  dydyydZZ  dyydZZ  dydyydZZ 

=  1 


Thus  we  have 


df d  d  +d2d  d  +d2d  d  -d2d2  -  d2d2  -  d2d2 

Uj'rUjyyUjzz  i  UjnUjxxUj  zz  >  ^  ^lAjxxUjyy  ^x^yz  UjyUjxz  ^z^xy 

j2  72  72 


^x^yy^zz  '  -y 
dxxdyy  T  dxxdzz  +  dyyd 


^xx^yy 

j2 


xx^yy  i  ^xx^zz  i  ^yy^zz  '  dXy  dxz  dyZ  dydyydzz  dzdyydzz 


=Gn 


— dzd  d  —  dzd  d  —  dzd  d  —  dzd  d 

UjxUjxxUjzz  UjzUjxxUjzz  UjyUjxxUjyy  UjxUjxxUjyy 


72  72 


+dzdz  Jrd2d2  4-  dzdz  4-  dzdz  +dzdz  +dzdz 

-TtLyLLyz  T  ^ztiyz  ~T  ^  x  ^  x  z  T-  ^z^xz  I  ^x^xy  '  ^y^xy 


& 

j2  72 


j2 


(7) 


Using  (7)  and  rearranging  the  rest  of  the  terms  in  (6)  we  obtain  G  =  Gv.  □ 
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Proposition  2.  Consider  a  C2  compact  surface  TcK"  (n  =  2,3)  of  codimension  1 
and  let  d  be  defined  as  in  (4) .  Define  the  closest  point  projection  map  Pp  os  in  (3)  for 
x  G  Rn.  For  \r)\  sufficiently  close  to  zero,  let  T^  be  the  rj  level  set  of  d 

rv  :=  {x  :  d(x)  =  r]}  .  (8) 


Define  the  Jacobian  Jv  as 

j  f  1  +  VKr)  ifn  =  2, 

^  (  1+2 r/H^  +  rfG^  if  n  =  3, 

where  is  the  signed  curvature  of  T^  in  2D,  and  Hv  and  Gv  are  its  Mean  curvature 
and  Gaussian  curvature  respectively  in  3D. 

Then  if  is  the  Jacobian  matrix  of  Pp  we  have 


r  01,  n  =  2, 

\  0102,  n  =  3, 


(9) 


where  cq,  cq  are  the  first  two  singular  values  of  the  Jacobian  matrix  . 

Proof.  The  distance  function  d  satisfies  the  property  d(x)  =  0  for  x  G  T.  Also,  since  T 
is  C2,  its  distance  function  d  belongs  to  (72(Rn,R);  see  e.g.  [1,  4].  It  follows  that  the 
order  of  the  mixed  partial  derivatives  does  not  matter.  In  addition,  the  normals  to  a 
smooth  interface  do  not  focus  right  away  so  that  the  distance  function  is  smooth  in 
a  tubular  neighborhood  T  around  T,  and  is  linear  with  slope  one  along  the  normals. 
Therefore  we  have 

|Vd|  =  1  for  all  x  G  T.  (10) 

The  third  important  fact  is  that  the  Laplacian  of  d  at  a  point  x  gives  (up  to  a  constant 
related  to  the  dimension)  the  mean  curvature  of  the  isosurface  of  d  passing  through 
x,  namely 

Ad(x)  =  (1  —  n)H(x),  (11) 

where  P(x)  is  the  Mean  curvature  of  the  level  set  {y  :  d( y)  =  d(x)}.  Differentiating 
(10)  with  respect  to  each  variable  gives  the  following  equations  in  three  dimensions: 

h'X h'XX  dydXy  +  dzdxz  0,  (12) 

dxdyx  ~\~  dydyy  -\-  dzdyZ  0,  (13) 

dxdZx  dydZy  -\-  dzdzz  0.  (14) 


In  particular  the  two  dimensional  case  can  be  derived  by  assuming  that  the  distance 
function  is  constant  in  z. 
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Two  dimensions.  In  that  case  the  Jacobian  matrix  dA  of  the  closest  point  pro¬ 
jection  map  is 

9Pr  f  1  dx  ddxx  ( dydx  T  ddyx ) 

( d'X  dy  T  ddXy ) 


dx 


1  ~  dy  —  ddyy 


Since  Schwartz’  Theorem  holds,  we  have  dxy  =  dyx  making  did.  a  real  symmetric 
matrix.  It  is  therefore  diagonalizable  with  eigenvalues  0  and  1  —  dAd.  Indeed,  we  have 


/  dx  (  1  dx  dy  )  d(  dxdxx  T  dydyx  )  \ 


dPr 

dx 


Vd  = 


=o  by  (10)  in  2D 

dy(  1-dl-dl  ) 


=0  by  (12)  in  2D 

d(y  dydyy  ~\~  dXdXy  ) 


=  0, 


\  =o  by  (10)  in  2D  =o  by  (13)  in  2D  J 


and  for  v  = 


dy 

dx 

dPr 

<9x 


V  = 


d/y  I  dydx  I  diy ddxx  dxdy  ddiXdiXy 

dn,dx  dyddxx  d  dx  dxddn 


dy 

dT 


-\-  d 


y^^xx 

dy  dxx  —  dxd. 


x^xy 


dydXy  dXdyy 


—  v  +  d 


(  dy^j  ( dydyy  ~\~  dXdXy^j  \ 

=0  by  (13)  in  2D 
A.d(dx^j  +  dxdxx  -\-  dydXy 


=o  by  (12)  in  2D  / 


=  (1  —  dAd)v. 


Since  ||v||  =  1,  v  is  an  eigenvector  corresponding  to  the  eigenvalue  A  =  1  — dAd.  Thus, 
for  x  such  that  d(x)  =  r/  we  have  that  the  eigenvalue  A  of  'dA-  satisfies 

A  =  1  —  r/Ad  =  1  +  r/Ky 

by  (11).  Since  1  +  r/Kv  >  0,  it  follows  that  A  coincides  with  the  singular  value  of  dA 
and  hence 

ai  —  l  +  riKv. 


Three  dimensions.  Since  for  \rj\  sufficiently  close  to  0  the  distance  function  is  C2 , 
the  Jacobian  matrix 


dPr 

dx 


i  dx  ddxx 
( dxdy  T  ddxy ) 
( dxdz  +  ddxz ) 


(dydX  "f  ddyp ) 
1  ~  dy  —  ddyy 

( dydz  ■*;  ddyz ) 


—  (dzdx  +  ddzx ) 
( dzdy  ddzy ) 

1  —  d2z  —  ddzz 


i 
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is  a  real  symmetric  matrix  which  is  diagonalizable  with  one  zero  eigenvalue  and  two 
other  eigenvalues  Ai  and  A2.  Indeed  using  (12), (13), (14)  and  (10)  we  can  show  that 


dPr 

dx 


Vd  =  0. 


Now  consider  x  such  that  d(x)  =  q.  Then,  the  characteristic  polynomial  ;\(A)  of 
is 

x(A)  =  —A  (A2  —  (2  —  r/Ad)X  —  Q) , 

where  Q  =  — Gr}r: f  +  qAd  —  1  with  Gn  dehned  in  (5).  Since  the  other  two  eigenvalues 
of  are  the  solutions  of  the  quadratic  equation  A2  —  (2  —  qAd)X  —  Q  =  0,  it 
follows  that 

A1A2  =  — Q  =  1  —  qAd  +  q2Grj  =  1  +  2qHv  +  q2G^. 

Since  1  +  2 qHv  +  q2Gv  >  0,  it  follows  that 


cri<72  =  1  +  2  qHv  +  q2Gv , 

where  cq  and  a2  are  singular  values  of  □ 

This  leads  to  the  following  proposition: 

Theorem  3.  Consider  T  a  curve  in  2D  or  surface  in  3D  with  C 2  boundaries  if  it  is 
not  closed,  and  define  d  :  Mn  1— »  M+  U  {0}  (n  =  2,  3)  to  be  the  distance  function  to  T 
with  Pr  :  Mn  1— >  T  the  closest  point  mapping  to  T.  Then 


f  u(x)dx  =  f  u(Pr(x)5e(d(x))E(x)dx, 

Jr  J  Rn 

where  Se  is  an  averaging  kernel  and  E  (x)is  defined  as 


E(x)  = 


di(x),  n  =  2, 

<7i(x)<72(x),  n  =  3, 


(15) 


where  crj(x)  ,  j  =  1,2,  is  the  j-th  singular  value  of  the  Jacobian  matrix,  evaluated 
at  x. 

Proof.  If  T  is  closed  we  combine  Equation  (2)  with  the  result  ./ (x)  =  E(x)  from 
Equation  (9)  of  Proposition  2. 

If  T  is  open  there  is  a  little  more  to  show  since  Equation  (2)  was  only  derived  for 
closed  manifolds.  Before  we  state  the  result,  it  is  necessary  to  understand  how  Tr) 
defined  in  (8)  (an  q— level  set  of  d)  looks  like  for  an  open  curve  in  two  dimensions  and 
for  a  surface  with  boundaries  in  three  dimensions. 

In  two  dimensions,  T,,  consists  of  a  flat  tubular  part  on  either  side  of  the  curve  and 
two  semi  circles  at  the  two  ends  of  the  curve.  See  Figure  1. 

In  three  dimensions  T  is  in  general  made  up  of  three  distinct  parts:  the  interior  part, 
the  edges  of  the  boundary  and  the  corners.  If  we  assume  that  T  has  N  edges  then 
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Figure  1:  An  example  of  an  open  curve  T  (black  curve)  and  its  77- level  set  r.,;  (red 
curve).  r„  consists  of  a  tubular  part  and  two  semi  circles  at  the  two  ends. 

we  can  write  T  =  r°  U  (U fLiEi)  U  (U f=1Ci),  where  T0  is  the  interior  of  T,  Ei  is  the 
i-th  edge  of  the  boundary  of  T  and  Ci  is  its  i  -th  corner.  In  that  setting  we  can  write 
r„  =  7„U  {uf=lT^)  U  (U^S^),  where  In  is  the  inside  portion  of  r^,  T-1  is  the  cylindrical 
part  of  r„  representing  the  set  of  points  located  at  a  distance  77  from  the  i  -th  edge 
Ei,  and  finally  S]1  is  the  spherical  part  of  F?/  representing  the  set  of  points  located  at 
a  distance  77  from  the  i  -th  corner  Ct .  See  Figure  2. 


1() 


Figure  2:  An  example  of  a  surface  with  boundaries  viewed  from  different  angles  (2(a), 
2(b)  and  2(c))  and  its  corresponding  77-level  set  r,;  viewed  from  the  same  angles  (2(d), 
2(e)  and  2(f)).  Figure  2(f)  shows  the  surface  and  F?; . 

In  both  cases  we  need  to  integrate  over  Yn  and  then  subtract  the  two  semi  circles  at 
the  two  end  points  of  the  curve  (in  two  dimensions)  or  subtract  the  portions  of  sphere 
at  the  corners  of  the  surface  and  the  portions  of  cylinders  at  the  edges  of  the  surface 
(in  three  dimensions).  However,  it  turns  out  that  the  subtraction  is  unnecessary  since 
E(x)  =  0  on  each  of  the  subtracted  pieces  as  shown  below. 

•  Two  dimensions.  On  the  semi-circle  around  the  end  point  of  a  curve,  the 
closest  point  mapping  is  constant  since  all  points  on  the  semi-circle  map  to 
the  end  point.  As  a  result,  the  singular  values  of  the  Jacobian  matrix  of  the 
closest  point  mapping  are  all  zeros  and  thus  E(x)  =  0  on  the  semi-circles  around 
the  end  points  of  a  curve. 
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•  Three  dimensions.  As  in  two  dimensions,  on  the  portions  of  sphere  around  a 
corner  point  of  a  surface,  the  closest  point  mapping  is  constant  and  thus  £(x)  = 
0.  On  the  portion  of  cylinders,  the  closest  point  mapping  is  constant  along  the 
radial  dimension  (one  of  the  principal  directions  or  singular  vector)  resulting  of 
the  singular  value  along  that  direction  to  be  zero.  Since  £(x)  is  the  product  of 
the  singular  values,  it  follows  that  £(x)  =  0  on  the  portion  of  cylinders  as  well. 

Consequently,  Equation  (15)  holds  for  any  C 2  curve  or  surface  with  C2  boundaries  of 
codimension  1.  □ 

2.2  Codimension  2 

We  consider  a  C2  curve  in  R3  denoted  by  T  and  let  y(s)  be  a  parameterization  by 
arclength  of  T.  We  denote  by  d  :  M3  R+  U  {0}  the  distance  function  to  T  and  let 
Pr  :  R3  |— t  r  be  the  closest  point  mapping  to  T.  We  consider  a  parameterization  of 
the  tubular  part  of  the  level  surface  for  rj  €  [0,  e]  defined  as 

x(s,0, 77)  :  7(s)  +  ?]cos#N(s)  +  r]sin6lB(s), 

where  T  =  ^,  N  and  B  constitute  the  Frenet  frame  for  7  as  illustrated  in  Figure  3. 
As  in  the  previous  section,  if  T  is  closed  then  d  is  the  signed  distance  function  to  T  . 


Figure  3:  Three  dimensional  curve  with  its  77-level  surface  F.,?  in  green  and  the  Frenet 
frame  at  a  point  on  rr). 

If  we  project  a  point  x  on  the  tubular  part  of  the  level  surface  r,;  defined  in  (8), 
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we  have  Pp(x(s,  0 ,  r/))  =  7(s).  If  L  is  the  length  of  the  curve  it  follows  that 

p2tt  pL  p2tt  pL 

/  /  g(Pr(x(s,  7)))|xs  x  x.g\dsd9  =  /  /  g(7(s))p(l  -  t)k(s)  cos  6)dsd6, 

J  0  J  0  JO 


0  J  0 

pL  p2tt 

=  rl  g(j(s))  /  (1  —  rjK cos  d)d6ds, 


(16) 


J  0  Jo 

=  27T7  J  g(j(s))ds. 

Note  that  the  tubular  part  of  the  level  surface  Fr)  does  not  contain  the  two  hemi¬ 
spheres  of  r„  which  are  located  at  the  two  end  points  of  the  curve  F.  Thus, 

[  g(Pr(x))dSx  =  2ivg  [  gds,  (17) 

Jrv\{Ciuc2}  Jr 

where  C\  and  Cb  are  the  two  hemispheres  of  the  level  surface  T,?  located  at  the  two 
end  points  of  the  curve  T.  Consequently,  for  sufficiently  small  e  and  by  the  coarea 
formula  we  obtain 


L 


9{lf{s))ds  =  2. 

l 

2n 


lo  ^7  Ir 

[  g(P f(x)) 

Jr  3 


r„\{Ciuc2} 

Kjd) 

d 


g{P f(x))  Ke(r))dg, 


X(c1uc2)-(x)dx, 


where  Ke  is  a  Cl  averaging  kernel  supported  in  [0,  e]  and  X(Ciuc2)c(x)  is  the  charac¬ 
teristic  function  of  the  set  ( C\  U  Cb)c.  In  our  numerical  simulations  we  consider  the 
kernel 


Kl’\v)  =  \  (i  -  cos  X[o,e](v)-  (18) 

Since  the  formulation  above  does  not  use  the  two  hemispheres  located  at  both  end 
points  of  the  curve,  in  order  to  integrate  over  the  tubular  part  of  rr)  only,  it  is  necessary 
to  subtract  the  integration  over  each  of  the  hemispheres  C\  and  Cb  .  The  result  can 
be  summarized  in  the  following  proposition: 

Proposition  4.  Consider  a  single  C2  curve  T  in  M3  parameterized  by  7 (s)  where  s  is 
the  arclength  parameter,  and  let  d  be  the  distance  function  to  T.  We  define  Ke  to  be 
a  C1  averaging  kernel  compactly  supported  in  [0,  e]  and  Pr  :  R3  1— >  T  to  be  the  closest 
point  mapping  to  T. 

If  g  is  a  continuous  function  defined  on  T  then  for  sufficiently  small  e  >  0  we  have 

f  9(l(s))ds  =  2-  [  ff(Pr(x))^^*^dx-  2  /  g(xv)r]Ke(r])dr /,  (19) 

Jr  27T  jK3  d(xj  J o 

where  x.v  is  a  point  on  a  sphere  of  radius  r/. 


10 


Note  that  for  the  computation  of  the  length  of  a  curve,  the  correction  terms  given 
by  integrating  over  both  C\  and  C'2  is 

f,KM  |S:N=  fW4,A  =  2,e. 

Jo  V  Jo  V 

This  simple  correction  is,  however,  not  suitable  for  more  general  cases  that  contain 
multiple  curve  segments  and  several  integrands.  We  shall  derive  a  more  elegant  and 
seamless  way  to  perform  such  correction  in  the  following  section. 

Now  if  we  consider  a  C2  curve  in  three  dimensions  and  let  Py  be  its  closest  point 
mapping,  we  have  the  following  proposition: 

Theorem  5.  Let  a(x )  be  the  nonzero  singular  value  of  and  let  g  be  a  continu¬ 
ous  function  defined  on  T.  If  y(s)  is  the  arclength  parameterization  of  T  and  if  e  is 
sufficiently  small  we  have 

\\  Ke(d)  , 
x))  -  cr(x)dx, 

where  d  is  the  distance  function  to  T. 

Proof.  Let  x  e  M3.  Then  there  exists  g  >  0,  such  that  xeT,. 

Case  1:  x  is  on  the  spherical  part  of  r,;  corresponding  to  the  ^-distance  to  either  of 
the  two  end  points  of  the  curve  T.  WLOG  we  assume  that  x  is  at  a  distance  g  from  the 
first  end  point  C\  parameterized  by  7(0).  The  result  is  the  same  if  x  is  on  the  other 
sphere,  i.e.  at  a  distance  g  from  the  other  end  point  C'2.  In  that  case,  Pr(x)  =  7(0) 
for  all  x  on  the  spherical  part  so  that  the  Jacobian  matrix  dPjJ^  =  0  .  Therefore,  for 
x  on  the  spherical  part  of  T^,  all  singular  values  of  the  Jacobian  matrix  are  zero. 

Case  2:  x  is  on  the  tubular  part  of  T^.  In  that  case,  if  we  use  the  Frenet  frame 
centered  at  the  point  x  =  x(s,  9,  rf)  €  T^  ,  we  can  write  x  in  the  new  coordinate 
system  (T,  N,  B)  as 


x  =  7(5)  +  uN  +  icB,  (21) 

where  u  =  0  is  the  coordinate  of  x  along  T,  v  is  the  coordinate  along  N  and  w  is  the 
coordinate  along  B.  Since  the  projection  Tr(x)  =  7(s)  does  not  depend  on  v  nor  w 
(since  the  plane  (N,  B)  is  normal  to  the  curve  T)  is  follows  that 

d-Pr(x)  _  <9Pr(x)  _ 
dv  dw 

On  the  other  hand,  we  have 

<9Pr(x)  d'y(s)  ds  d'y(s)  ds^ 

du  du  du  ds  du  ’ 
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where  ||  is  the  variation  of  the  arclength  parameter  s  with  respect  to  u  when  the  point 
x  is  moving  on  Fr/  along  the  tangential  direction  T.  Since  u  is  the  arclength  parameter 
along  the  tangential  direction  T,  it  follows  that  we  have  a  unit  speed  parameterization 
along  T  giving  the  identity 


In  addition, 


dx  d^(s)  N  B 

&  =  jr  +> + -y*  __ 

=  T  —  kvT  +  tv  B  —  ™N 
=  (1  —  kv)  T  —  twN  +  tv B, 

where  k  is  the  curvature  of  F  at  y(s)  and  r  is  the  torsion  of  the  curve  T  at  the  point 
7(s).  Since  the  level  surface  r,;  is  a  tube  of  radius  //,  its  intersection  with  the  normal 
plane  (N,  B)  is  a  circle  of  radius  r/.  Hence  if  we  use  polar  coordinates  on  the  normal 
plane,  we  obtain  v  =  r/  cos  6  and  w  =  r/  sin  0  .  It  follows  that 


KTj  COS  9. 


Consequently  we  have 


and 


<9x  ds  dx 

- T  =  1  = - 

du  du  ds 


ds 

du 


ds 

—  (1-K77COS0), 


1  —  UT]  COS  6 

Therefore,  in  the  Frenet  frame,  the  Jacobian  matrix  of  the  closest  point  projection 
map  can  be  written  as 


i 


dPT 

dx 


l  —  KTj  COS  6 

0 


0 


where 

Since 


l  —  KTJ  COS  6 


1  —  KT]  COS  t 


1  —  KTl 

Theref* 


is  the  nonzero  eigenvalue  of  the  Jacobian  of  the  closest  point  mapping. 
>  0,  it  is  also  its  singular  value  a(x). 


'ore  we  have 


a(x) 


1 


1— k;?7Cos  6 


if  x  is  on  the  spherical  part  of 
if  x  is  on  the  tubular  part  of  Yrj . 


(22) 
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Now  using  (16)  and  (17)  we  obtain 


g(Pr(x))er(x)dSx  = 


g(Pr(x))<r(x)dSx 


'rv\{Ci\jc2} 

nL 

5f(Pr(x))c7(x)|xs  X  Xg\dsd9 

= r  fLgh  (s))d-nK{s>ms%sdo 

Jo  Jo 

=  2tt?7  /  g(^(s))ds 

Jo 


1  —  t)k(s)  cos  9 


It  follows  that  for  Ke  a  C1  averaging  kernel  compactly  supported  in  [0,  e],  for  sufficiently 
small  e  and  by  the  coarea  formula,  we  have 


J-  f  ~  f  g{Pri?0)<y{?0Ke{ri)dg 

jo  /  j r ^ 

h  L3(PT{x))Eldllj{x)dx' 


□ 


3  Numerical  simulations 

In  this  section  we  investigate  the  convergence  of  our  numerical  integration  using  simple 
Riemann  sum  over  uniform  Cartesian  grids.  In  our  computations  we  use  the  cosine 
kernel 

Ke°S(v)  =  X[0,e](h)^  (l  +  COS 

for  integration  on  surfaces  of  codimension  1,  and  the  kernel  K^1  defined  in  (18)  for 
codimension  2.  Unless  stated  otherwise,  the  singular  values  are  computed  from  the 
matrix  whose  elements  are  computed  by  the  standard  central  difference  approximations 
of  the  Jacobian  matrix  With  these  compactly  supported  kernels,  formulas  (15) 
and  (20)  can  be  considered  integration  of  functions  defined  on  suitable  hypercubes, 
periodically  extended.  In  such  settings,  simple  Riemann  sums  on  Cartesian  grids  are 
equivalent  to  sums  suing  Trapezoidal  rule,  and  the  order  of  accuracy  will  be  related  in 
general  to  the  smoothness  of  the  integrands;  exception  can  be  found  when  the  normals 
of  the  surfaces  are  rationally  dependent  on  the  step  sizes  used  in  the  Cartesian  grids. 

3.1  Integration  of  codimension  one  surfaces 

We  tested  our  numerical  integration  on  two  different  portions  of  circle,  a  torus,  a 
quarter  sphere  and  a  three  quarter  sphere.  We  computed  their  respective  lengths  or 
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surface  areas  by  integrating  the  constant  1  over  the  curve  or  surface.  Each  of  these 
tests  were  designed  to  exhibit  the  convergence  rate  of  our  formulations  on  cases  with 
varying  difficulty.  In  particular,  the  convergence  rate  of  our  formulation  depends  on 
the  smoothness  of  the  closest  point  mapping  inside  the  tubular  neighborhood  of  the 
surfaces. 

The  results  for  the  portions  of  circle  are  given  in  Tables  1  and  2.  In  the  first 
convergence  studies  (Table  1)  the  line  where  the  closest  point  mapping  has  a  jump 
discontinuity  is  parallel  to  the  grid  lines.  In  this  case  we  see  a  second  order  convergence 
rate  using  central  differencing  to  compute  the  Jacobian  matrix  In  the  second  test 
case  however,  the  portion  of  circle  is  chosen  so  that  the  line  where  the  closest  point 
mapping  has  a  jump  discontinuity  is  not  parallel  to  the  grid  lines.  In  that  case  the 
normal  to  the  curve  is  rationally  dependent  on  the  step  size  of  the  Cartesian  grid  and 
the  convergence  rate  reduces  to  first  order  even  though  we  used  central  differencing  to 
compute  We  note  that  in  these  two  tests,  we  chose  e  (the  half  width  of  the  tubular 
neighborhood  around  the  curve)  small  enough  so  that  the  line  where  the  closest  point 
mapping  is  discontinuous  is  outside  of  it. 

In  three  dimensions  we  first  tested  our  method  on  a  torus  (closed  smooth  surface) . 
The  results  for  the  torus  are  reported  in  Table  3.  In  this  case  the  closest  point  mapping 
is  very  smooth  and  we  see  third  order  convergence  when  using  a  third  order  difference 
scheme  to  approximate  For  surfaces  with  boundaries  we  tested  the  method  on  a 
quarter  sphere  and  a  three  quarter  sphere.  The  three  quarter  sphere  case  is  illustrated 
in  Figure  4. 


Figure  4:  The  three  quarter  sphere  (4(a))  and  its  corresponding  //-level  set  rr/  (4(b)). 

The  reason  for  choosing  these  two  cases  is  because  the  closest  point  mapping  has 
a  different  degree  of  smoothness  for  each  of  these  surfaces.  For  the  quarter  sphere  the 
closest  point  mapping  is  smooth  enough,  but  for  the  three  quarter  sphere,  the  tubular 
neighborhood  around  the  surface  contains  the  line  where  the  closest  point  mapping  has 
a  jump  discontinuity.  In  that  latter  case,  it  is  therefore  necessary  to  use  an  adequate 
one  sided  discretization  to  compute  accurately.  The  one-sided  discretization  that 
we  used  is  reported  in  Section  3.3.  The  test  for  the  quarter  sphere  still  uses  central 
differencing  to  compute  The  results  for  the  portions  of  sphere  are  reported  in 

Tables  4  and  5. 
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Table  1:  Relative  errors  in  the  numerical  approximation  of  the  length  of  a  planar 
curve,  which  is  a  portion  of  circle  of  radius  R  =  0.75  centered  at  0.  The  width  for  the 
tubular  neighborhood  of  the  curve  is  e  =  0.2.  In  this  computation,  the  closest  point 
mapping  has  a  jump  discontinuity  along  a  straight-line  which  is  arranged  to  be  parallel 
to  the  grid  lines. 


n 

Relative  Error 

Order 

64 

2.7994  x  10"4 

- 

128 

7.0665  x  10"5 

1.99 

256 

1.7187  x  10"5 

2.04 

512 

4.2719  x  10“6 

2.01 

1024 

1.0636  x  10“6 

2.01 

2048 

2.6567  x  10“7 

2.00 

4096 

6.6045  x  10“8 

2.01 

8192 

1.6513  x  10"8 

2.00 

Table  2:  Relative  errors  in  the  numerical  approximation  of  the  length  of  a  planar  curve, 
which  is  a  portion  of  circle  of  radius  R  =  0.75  centered  at  0.  The  width  for  the  tubular 
neighborhood  of  the  curve  is  e  =  0.2.  In  this  computation,  the  jump  discontinuity  of 
the  closest  point  mapping  is  not  parallel  to  the  grid  lines. 


n 

Relative  Error 

Order 

64 

3.7159  x  10“5 

- 

128 

2.5786  x  10“7 

7.17 

256 

4.2361  x  10“6 

-4.04 

512 

3.2246  x  10"6 

0.39 

1024 

1.8876  x  10"6 

0.77 

2048 

1.0132  x  10-7 

0.90 

4096 

5.2372  x  10-7 

0.95 

8192 

2.6615  x  10“7 

0.98 
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Table  3:  Relative  errors  in  the  numerical  approximation  of  the  surface  area  of  a  torus 
centered  at  0.  The  distance  from  the  center  to  the  tube  that  form  the  torus  is  R  =  0.75 
and  the  radius  of  the  tube  is  r  =  0.25.  In  this  computation,  we  summed  up  grid 
points  that  are  within  e  =  0.2  distance  from  the  surface.  The  Jacobian  matrix  is 
approximated  by  a  standard  third  order  accurate  differencing. 


n 

Relative  Error 

Order 

32 

6.2030  x  10“3 

— 

64 

1.8073  x  10“4 

5.10 

128 

6.6838  x  10"6 

4.76 

256 

4.1530  x  10-7 

4.01 

512 

5.0379  x  10"8 

3.04 

Table  4:  Relative  errors  in  the  numerical  approximation  of  the  surface  area  of  a  quarter 
sphere  with  radius  R  =  0.75  centered  at  0.  In  this  computation,  we  summed  up  grid 
points  that  are  within  e  =  0.2  distance  from  the  surface.  We  used  the  standard  central 
difference  scheme  to  compute  each  entry  of  the  Jacobian  matrix  ^p1 . 


n 

Relative  Error 

Order 

32 

9.2825  x  10“3 

— 

64 

1.8365  x  10"3 

2.34 

128 

2.7726  x  10"4 

2.73 

256 

7.1886  x  10"5 

1.95 

512 

1.4811  x  HT5 

2.30 

Table  5:  Relative  errors  in  the  numerical  approximation  of  the  surface  area  of  a  3/4 

sphere  with  radius  R  =  0.75  centered  at  0  (this  is  the  portion  of  a  sphere  that  misses 

half  of  a  hemisphere).  In  this  computation,  we  summed  up  grid  points  that  are  within 

e  =  0.2  distance  from  the  surface.  Due  to  this  setup,  the  closest  point  mapping 

has  a  discontinuity  that  stems  out  from  the  boundary  of  the  surface.  We  used  the 

discretization  described  in  Section  3.3  to  compute  each  entry  of  the  Jacobian  matrix 
dPr 
dx  ■ 


n 

Relative  Error 

Order 

32 

1.1726  x  10"2 

— 

64 

1.1733  x  10"3 

3.32 

128 

9.1325  x  10“4 

0.36 

256 

3.8238  x  10“4 

1.26 

512 

7.8308  x  10“5 

2.29 
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3.2  Integrating  along  curves  in  three  dimensions 

In  codimension  2,  we  tested  our  numerical  integration  on  a  coil  wrapped  around  the 
helix  defined  parametrically  as 

x(i)  =  (r  cos(t),r  sm(t),bt) , 

with  r  =  0.75  and  b  =  0.25.  The  coil  is  then  wrapped  around  the  helix  at  a  distance  of 
0.2  from  the  helix.  As  our  test  case,  we  computed  the  length  of  the  coil  by  integrating 
1  along  the  curve.  The  results  are  reported  in  Table  6. 


Figure  5:  The  coil  and  one  of  the  level  sets  of  the  distance  function  to  the  coil  used  in 
the  reported  numerical  simulations. 


Table  6:  Relative  errors  in  the  numerical  approximation  of  a  coil  wrapped  around  a 
helix.  In  this  computation,  we  used  a  constant  width  for  the  tubular  neighborhood 
e  =  0.1  and  took  the  averaging  kernels  to  be  AT4,1  defined  in  (18). 


n 

Relative  Error 

Order 

60 

5.5078  x  10"3 

— 

120 

1.1476  x  10"3 

2.63 

240 

2.3409  x  10"4 

2.29 

480 

3.7166  x  10“5 

2.66 

3.3  One-sided  discretization  of  the  Jacobian  matrix 

Here  for  completeness,  we  describe  the  one-sided  discretization  used  in  computing 
results  reported  in  Table  5. 

For  simplicity,  we  will  describe  the  one-sided  discretization  for  a  uniform  Cartesian 
grid,  namely  for  Pr(xjj)  =  with  Xjj  =  ( ih,jh ),  i,j  G  Z  and  h  >  0  being 
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the  step  size.  The  Jacobian  matrix  will  be  approximated  by  simple  finite  differences 
defined  below: 


( UX)i,j  ( Uy)iJ  \ 

(Uk i  (Vv)i,j  ) 


The  discretization  of  U  and  V  have  to  be  defined  together  because  the  two  functions 
are  not  independent  of  each  other.  With 


(Ux)ij  ±2^  Ui±2,j)  5 

and  the  smoothness  indicator 


we  define 


Sy  =  ^(tfy)  :=  A+A ~Ui±hj 

(Ux)ij  := 


f (^+)«.  ^  is;i  <  15--1, 

1  (U~)itj,  otherwise, 
and  iyx)i,j  is  defined  according  to  the  choice  of  stencil  based  on  S±(Uhj) 


lvv..=  J(U+)«.  if|SJI<|Skl. 

otherwise. 

The  discretization  of  U„  and  V,,  is  defined  similarly  with  the  choice  of  the  stencil  based 
on  SHV,j). 


4  Summary 

In  this  paper,  we  presented  a  new  approach  for  computing  integrals  along  curves  and 
surfaces  that  are  defined  either  implicitly  by  the  distance  function  to  these  manifolds 
or  by  the  closest  point  mappings.  We  are  motivated  by  the  abundance  of  discrete 
point  sets  sampled  from  surfaces  using  devices  such  as  LIDAR,  the  need  to  compute 
functionals  defined  over  the  underlying  surfaces,  as  well  as  many  applications  involving 
the  level  set  method  or  the  use  of  closest  point  methods. 

Contrary  to  most  other  existing  approximations  using  either  smeared  out  Dirac 
delta  functions  or  locally  obtained  parameterized  patches,  we  derive  a  volume  integral 
in  the  embedding  Euclidean  space  which  is  equivalent  to  the  desired  surface  or  line 
integrals.  This  allows  for  easy  construction  of  higher  order  numerical  approximations 
of  these  integrals.  The  key  components  of  this  new  approach  include  the  use  of  singular 
values  of  the  Jacobian  matrix  of  the  closest  point  mapping,  which  can  be  computed 
easily  to  high  order  even  by  simple  finite  differences. 
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